Symmetry breaking and coarsening of clusters in a prototypical driven granular gas 
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Granular hydrodynamics predicts symmetry-breaking instability in a two-dimensional (2D) en- 
semble of nearly elastically colliding smooth hard spheres driven, at zero gravity, by a rapidly 
vibrating sidewall. Super- and subcritical symmetry-breaking bifurcations of the simple clustered 
state are identified, and the supercritical bifurcation curve is computed. The cluster dynamics pro- 
ceed as a coarsening process mediated by the gas phase. Far above the bifurcation point the final 
steady state, selected by coarsening, represents a single strongly localized densely packed 2D cluster. 

PACS numbers: 45.70.Qj 



Introduction. Rapid granular flow continues to at- 
tract a great deal of attention of physicists Q. Among 
most fascinating phenomena here is clustering: nuclc- 
ation and growth of dense granular clusters, surrounded 
by dilute granular gas in "freely-cooling" and driven 
H [| [| ||, 0] granular gases. Clustering can be viewed 
as a variant of thermal condensation instability, encoun- 
tered also in gases and plasmas that cool by their own 
radiation H. This analogy brings about the question of 
universality of structure formation in ensemble of par- 
ticles with energy losses of different nature. A related, 
largely unexplored issue is cluster coarsening. Reviewing 
different cluster-forming granular systems with a fixed 
number of particles J|, 0, || , one notices that coarsening is 
ubiquitous. In most of these systems a single cluster usu- 
ally survives after transients die out. The present work 
addresses novel issues of pattern formation and coarsen- 
ing in cluster-forming driven granular gases. We shall 
consider a simple (indeed, prototypical) system: a 2D 
ensemble of inelastically colliding smooth hard spheres, 
confined in a rectangular box and driven by a rapidly vi- 
brating side wall at zero gravity. The other three walls 
are assumed elastic. Though this and related systems 
were investigated theoretically (|, ||, |ll| and exper- 
imentally E3, it has been recognized only recently 
that they exhibit nontrivial pattern-forming properties 
p3| , [l4| |. These properties will be in the focus of this 
work. We identify the character of symmetry-breaking 
bifurcations of the simple clustered state and show that, 
depending on the control parameters, both super- and 
subcritical bifurcations can occur. We find that selection 
of the final steady state occurs via cluster coarsening dy- 
namics, mediated by the gas phase. Far above the bifur- 
cation point, only one densely packed 2D cluster survives. 
We shall conclude that cluster selection by coarsening is 
universal in cluster-forming granular flows with a fixed 
number of particles. 

Model. Let N 3> 1 identical smooth hard spheres of 
diameter d and mass m = 1 are rolling without friction 
on a smooth horizontal surface of a rectangular box with 
dimensions L x x L y . The number density of grains is 
n(x,y,t), the granular temperature is T(x,y,t). For a 



submonolayer coverage n < n c — 2/{-\fid 2 ), the (hexago- 
nal) close-packing density. Three of the walls are immo- 
bile, and grain collisions with them are assumed elastic. 
The fourth wall (located at x — L x ) is rapidly vibrat- 
ing, x = L x + Acosut, and supplies energy to the sys- 
tem. The energy is being lost via inelastic hard-core grain 
collisions characterized by a constant normal restitution 
coefficient r. Our crucial assumption will be a strong in- 



equality 1 



« 1. In the quasielastic limit, and for 



small Knudsen numbers, the Navier-Stokes granular hy- 
drodynamics is expected to be reasonably accurate, even 
for large granular densities, as long as the granulate is in 
a fluidized state [||. Therefore, the steady states of the 
system can be described by the coarse-grained momen- 
tum and energy balance equations 



p = const , V ■ (k VT) = / , 



(1) 



where p is the pressure, n is the thermal conductivity 
and / is the rate of energy losses. To make the full use of 
hydrodynamics, we shall work in the parameter regime 
when the energy supply from the vibrating wall can be 
represented as a hydrodynamic heat flux. This requires 
A <C I, where I is the mean free path of the particles at the 
vibrating wall. We shall also assume a double inequality 
T 1 / 2 /I < lu < T X I 2 IA. The left inequality guarantees 
the absence of correlations between successive collisions 
of particles with the vibrating wall. The right inequality 
makes the calculation simpler, but it is not crucial. The 
resulting energy flux is (ll| [l^] 



q = k OT/dx = (2/tt) 1 / 2 A 2 to 2 n T 1 ' 2 



(2) 



To close the hydrodynamic model, one needs constitutive 
relations (CRs) : p, k and I in terms of n and T. The 
CRs have been derived systematically only in the dilute 
limit. Reasonably accurate CRs in the whole range of 
densities can be obtained by employing free volume ar- 
guments close to the dense-packing limit, interpolating 
between the high- and low-density limits and finding the 
fitting constants by comparing the results with particle 
simulations [Q, ^6[. We shall use the CRs suggested by 
Grossman et al. B because of their relative simplicity. 
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A special investigation |Q showed that the stability di- 
agram does not change much if one uses instead the CRs 
derived by Jenkins and Richman [ jl7| . 

Eqs. (|l|) can be rewritten in terms of one variable: the 
(scaled) inverse density z(x,y) = n c /n(x,y) ||, We 
introduce scaled coordinates r/L x — > r so that the box 
dimensions become 1 x A, where A = L y /L x is the box 
aspect ratio. We have 



V • (F(z)'Vz) = £Q(z). 
Introducing tp — J Q F(z') dz' , we rewrite it as 

The boundary conditions are 

= 



dtp 




dtp 




dtp 




dx 


x=0 


dy 


y=0 


dy 


y=A 



(3) 



(4) 



(5) 



and (see Ref. Q) 



dtp 
dx 



CH[iP(l,y)] 



JjQ(tP)dxdy 
o o 

jH(tP(l,y))dy 
o 



(6) 



Here Q(tp) = Q [z(tp)} and H(ip) = H [z(ip)], while 
F, G, H and Q are functions of z only; they are given 
in Ref. |[3) . In the rest of the paper the symbol ~ will be 
omitted. Conservation of the total number of particles 
yields an equation for the area fraction / = N / (L x L y n c ): 



A 1 



1 

A 



dxdy 
^P) 



(J) 



The steady state problem is fully determined by three 
scaled parameters: C — (32/37) (L x /d) 2 (1 — r 2 ) (where 
7 ~ 2.26), the area fraction / and aspect ratio A. Notice 
that the steady-state density distributions are indepen- 
dent of A and uj. This is in contrast to the case of a non- 
zero gravity, where the gravity acceleration, combined 
with Au> 2 , yields a governing parameter. 

Strip state. The basic state of the system is the "strip 
state" : a laterally symmetric cluster located at the wall 
x = 0. The physics of the strip state is simple. Due to 
inelastic collisions, the granular temperature goes down 
with increasing distance from the driving wall. To main- 
tain the pressure balance, the granular density should in- 
crease with this distance, reaching the maximum at the 
opposite (elastic) wall. The strip state is described by 
the y-independent solution of Eqs. (^) and (@); we shall 
call it z = Z{x) that corresponds to tp = ^>(x). Notice 
that Eq. (j^) holds automatically in ID |l3). An example 
of the strip state is shown in Fig. 2 (bottom left). A 
similar state was observed in experiment |6| . 



Symmetry-breaking instability. The ID strip state 
gives way, by a symmetry-breaking bifurcation, to 2D 
clustered states. The bifurcation point can be found by 
linearizing Eqs. (^)-(0) around tp = ty(x). In the frame- 
work of time- dependent hydrodynamics, this corresponds 
to marginal stability of the strip state with respect to 
small perturbations along the strip. We write 



tp(x, y) = ~$f(x) + (fik{x) exp(ifcy) + c.c. 



(8) 



and, after linearization, arrive at a linear eigenvalue prob- 
lem for k — k c (f): 



Lpl - CQm, ip k - k 2 c ip k = . 

^(o) = o, 



(9) 



(10) 



/ Q(*(x))dx 



Here and below (. . .)*(x) = [F 1 d(. . .)/dz\ \ z _ Z ( x y anc ^ 
(...) stands for any function. Solving the eigenvalue 
problem for a given C, we obtain the marginal stabil- 
ity curve k = k c (f) and corresponding eigenfunctions 
fk{x). The modes with k < k c (f) are unstable. At 
fixed £, the instability occurs for f\{C) < f < /2OC), 
where fc c (/i) = &c(/2) = [|3|. The driving force of 
the instability is negative effective lateral compressibil- 
ity of the gas 0]. Figure 1 gives an example of the 
marginal stability curve in terms of the minimum aspect 
ratio A c (/) = ir/k c (f), at which the instability occurs. 
2 
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FIG. 1: The critical aspect ratio A for the instability and 
parameter B determining the bifurcation curve (J2TJ) versus 
the area fraction / for £ = 1.25 ■ 10 4 . The vertical dash-dot 
line corresponds to / = 0.0378. 

Bifurcation curve. To determine the nature of the bi- 
furcation (sub- or supercritical) and compute the bifur- 
cation curve, one should go to the second order of the 
perturbation theory. We can write 

tp(x,y) = $(x) + ^2 Vn{x) exp(inky) , (12) 

n 

where ip- n (x) = (p* n (x), and assume that cpo ~ ip\, 
(f2 *~ (pi, (pa ~ if 1, etc. Therefore, we need to take into 
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account only the terms n — 0, ±1 and ±2. This yields 
the following linear equations: 



C 



ip" — C Qsa <p\ — k 2 ipi = (k 2 - k 2 c ) ip+ 



1 2 

((pop + (p2<p*) + 2 Q*** <p m 



(13) 



(p'2 - £ (fi2 - 4,k 2 ip 2 = ^ & Qi"f V 2 > 
and the boundary conditions: 

^(0) = ^(0)=^ 2 (0) = 0, 



(14) 



(15) 



(16) 



Hqj 

IT 



<pi I Q dx 



J Q dx 



^—^ — f Hm-q, (tp (p + (f 2 v*) + <P \<f\ 2 



i l 



\ip\ dx ip 



Hq, 
H 2 



X 

(^H^ (p a + Hq,^ \ip\ 2 ^J jQdx 



Hq, Hq,q, 2 

^2 + -777T m 



H 



2H 



(17) 



>dx, (18) 



where ip = AY(x) is a properly normalized solution of 
Eqs. (j^)- (JTTJ) , see below. As the boundary condition for 
(po at x = 1 is fulfilled automatically, one more condition 
is needed. For a fixed /, this condition is supplied by Eq. 
(0): 



^■dx = 2 1 ( 1 



Z 2 F 



\Z 3 F 2 2Z 2 F 2 



y '"' * M 2 dx. (19) 



the solution of Eqs. (||) and ( ^Cj| ) obeying the normal- 
ization Y(0) = 1. This yields A (k 2 - k 2 ) = CA\A\ 2 , 
where C = const. The trivial solution A = describes 
the strip state, while the nontrivial one, k 2 — k 2 = C\A\ 2 
describes the bifurcated state. C > (< 0) corresponds 
to supercritical (subcritical) bifurcation. The solvability 
condition is a generalization of the standard "orthogonal- 
ity" condition, or the Fredholm alternative |ll| . It yields 
C explicitly in terms of definite integrals of solutions of 
the homogeneous forms of Eqs. (g), @ and @ that 
can be found numerically, see Appendix. We present here 
the resulting bifurcation curve for Y" c , the (normalized) 
y-coordinate of the center of mass of the granulate: 



Io dx I-A/2 v z ld y 



dy 



(20) 



where we shifted the y-coordinate: y + A/2 — > y. Let 
the aspect ratio of the system A be slightly larger than 
A c = 7r/fc c (/) so that only one mode, with k = 7r/A, is 
unstable. The bifurcation curve takes the form 



\YJ = 



1/2 



(21) 



where B = Cf 2 /(2k 2 Jf) and /i = 2 /„* YZ^F^dx. 
Eq. (|l]) assumes B > 0: a supercritical bifurcation. Fig- 
ure 1 shows B(f) for £ = 12, 500. We find that B > 
on an interval of / that lies within the instability interval 
(fi, f2)- Closer to the points /i and f 2 we obtain B < 
which indicates subcritical bifurcation. Subcritical bifur- 
cations close to the high-density instability border were 
previously observed by solving numerically the nonlinear 
steady state equations (13)- (0) f0|. 
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FIG. 2: Steady states found in time-dependent hydrodynamic 
simulations for £ = 1.25 • 10 4 , / = 0.0378 and A = 1.2 
(bottom left), 1.3 (top left) and 3.0 (right). The left wall of 
the box is the driving wall. 



The solvability condition for Eq. ([14|) yields the bifurca- 
tion curve: a relation between the amplitude of <p\ (we 
call it A) and k 2 — k 2 . One way to define A is the fol- 
lowing: (fi(x) = AY(x) + A\A\ 2 Sipi(x), where Y(x) is 



Time- dependent hydrodynamic simulations: bifurca- 
tions and coarsening. We performed a series of time- 
dependent hydrodynamic simulations, to verify the bi- 
furcation theory and to follow the cluster dynamics at 



4 



large aspect ratios. The full time-dependent hydrody- 
namic equations were solved with the same constitutive 
relations and boundary conditions as those used in our 
steady state analysis. Instead of the shear viscosity in the 
Navier-Stokes equation we accounted for a small model 
friction force — nv/r, where v is the hydrodynamic veloc- 
ity. An extended version of the compressible hydro code 



o o o 



o o o 



o o o 



VULCAN 1 19 was employed. 




FIG. 3: Bifurcation curve Y C (A) predicted by Eq. (pi (line) 
and found in hydrodynamic simulations (squares) for C = 
1.25 • 10 4 and / = 0.0378. 

We fixed C = 1.25 ■ 10 4 and / = 0.0378 and varied 
A. The initial scaled density included a zero mode cor- 
responding to the fixed / plus small-amplitude random 
noise. Figure 2 shows the final states for different aspect 
ratios A. For C and / used, the marginal stability theory 
predicts A c ~ 1.28 (see Fig. 1). Indeed, the strip state 
observed at A = 1.2 (Fig. 2, bottom left) gives way to 
a slightly asymmetric state at A = 1.3 (Fig. 2, top left). 
Far above the bifurcation point the final state represents 
a densely packed 2D cluster ("island"), located in a cor- 
ner (Fig. 2, right). This implies that all but one of the 
multiple 2D steady state solutions found earlier (chains of 
islands periodic in the y-direction) |l3| are unstable. The 
stable steady state selected by the coarsening dynamics 
is the one with the maximum possible period, equal to 
twice the lateral dimension of the system. 

Figure 3 shows the bifurcation curve Y C (A) predicted 
by Eq. (^l|), and measured in the simulations after tran- 
sients die out. Excellent agreement is obtained for not 
too large supercriticalities. Close to the bifurcation point 
we observed exponential slowdown as expected. 

Now we present the simulation results on the cluster 
dynamics and selection. For large A the dynamics involve 
two stages (see Fig. 4). During the first stage, several 
clusters nucleate at the wall opposite to the driving wall 
(Fig. 4A). Their number is of the order of A/A c which 
apparently corresponds to the maximum linear growth 
rate of the instability versus k. At the slower second stage 
the clusters become denser. As they compete for the 
material, their number goes down, and only one densely 
packed cluster ("island") finally survives, always in a cor- 
ner (Fig. 4e). The clusters interact mostly through the 
gas phase, similarly to Ostwald ripening in phase order- 
ing systems with conserved order parameter, controlled 



FIG. 4: Time history of the density field for C = 1.25 ■ 10 4 , 
/ = 0.0378 and A = 5 at scaled times 650 (A), 2, 100 (B), 
2,850 (C), 3,250 (D), and 7,000 (E). The left wall is the 
driving wall. Notice the change of color code with time. 



by gasdynamics g, |2(J. We also observed, for a different 
realization of noise in the initial conditions, direct coales- 
cence of transient clusters. However, the resulting single 
cluster is always the same in simulations with the same 
C, f and A. 

In summary, we used hydrodynamics to determine the 
character of symmetry-breaking bifurcations in a proto- 
typical driven granular gas and compute the supercrit- 
ical bifurcation curve. We found the selected steady 
state with broken translational symmetry and showed 
that selection is made via a coarsening process similar 
to Ostwald ripening. It appears that cluster selection by 
coarsening is a universal selection mechanism in cluster- 
forming granular flows with a fixed number of particles. 

This research was supported by the Israel Science 
Foundation, founded by the Israel Academy of Sciences 
and Humanities, and by the Russian Foundation for Ba- 
sic Research (grant No. 02-01-00734). 

Appendix. Solvability condition. Here we present 
the solvability condition for the boundary value prob- 
lem described by the linear equation (|l4j) and bound- 
ary conditions ( |l6| ) and (17). This solvability conditions 
yields the constant C that enters the bifurcation equation 
k 2 — k 2 — C ' \A\ 2 . Consider the following problem: 



w"{x)+P(x)w(x) = f(x). 

w'(0) = 0, 
w'(l)+aw(l) = g, 



(22) 
(23) 
(24) 



where g is a parameter; and assume that the homoge- 
neous variant of this problem, 



w"(x) + P(x)w(x) = 0, 
w'{0) = 0, 
w'(l) + aw(l) = 



(25) 
(26) 
(27) 
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has a nontrivial solution. The problem (E2|)-(B4|) will have 
no solution unless there is a relation between f(x) and 
g. What is the relation? To answer this question, wc 
construct the solution of the problem (22)-(E4fi using two 



fundamental solutions wq and w\ of the homogeneous 
equation (p5|). We choose these solutions as follows. Let 



wq obeys the boundary conditions (26) and 



iflo(O) = 1. 



(28) 



In view of the above assumption, Eq. ( |27| ) is fulfilled 
automatically. The second fundamental solution wi of 
Eq. ( p5| ) obeys the following boundary conditions (again 
at x = 0): 



wi(0) = 0, w[(0) = l. 



(29) 



Notice that each of the solutions wq and u>i is defined 
by an initial value problem, rather than by a boundary 
value problems. 

The general solution of Eqs. (p2|)-(p4]) can be writ- 
ten as Cq(x) wq(x) +Ci(x) w\(x), where functions Co (a:) 
and C\(x) are yet unknown, and we can demand 
C (x) wo(x) + C[ (x) w\ (x) = 0. Straightforward calcula- 
tions yield the solvability condition 



g = [w[(l)+aw 1 (l)] J fwo 



dx . 



(30) 



In the particular case of g — Eq. ( p0| ) is reduced to the 

l 

standard orthogonality condition J fwo dx = . 

o 

Applying the solvability condition (3C ) to Eq. (fhfl) with 
the boundary conditions ( fl6| ) and (17), we obtain the 
following relationship: 



(iPQip + tp 2 <p ) + OEJ v m 



H 



2H 



Q dx+ 



--JfV ( J Vo dx + J Q** \tp\ 2 dx 



H3, Hq/Hxj/^ 



H 2 



H 2 



l 

(p|(^| 2 ^ J ' Qdx 



x = l 



± i 
Y^-C^YnJ Qdx Jy w 



x=l 



(k 2 - k 2 ) V + 



t 1 2 

C I (tp <P + <Pz<f '*) + 2 Q*** V Ivl 



dx, (31) 



where Yq(x) is solution of Eq. (^|) with the boundary 
conditions 



F (0) = 1, F '(0) = 0, 



(32) 



and Y\ (x) is the solution of Eq. (|^) with the boundary 
conditions: 



yi(o) = o, y((o) = i. 



(33) 



Notice that each of the solutions Yq(x) and Y\{x) is de- 
fined by an initial value problem, rather than by a bound- 
ary value problems. This makes the numerical solution 
straightforward. To complete the calculation, we solve 
numerically the initial-value problems for the linear dif- 
ferential equations for ip, ipo and tp2 and compute the 
definite integrals entering Eqs. ([l9]) and (|l]). The final 
result is given by Eq. (pi]) . 
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